// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtIsovolumeFilter.C
// ************************************************************************* //

#include <avtIsovolumeFilter.h>

#ifdef HAVE_LIBVTKM
#include <avtVtkmDataSet.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/filter/contour/ClipWithField.h>
#include <vtkm/filter/field_conversion/PointAverage.h>
#endif

#include <vtkCellData.h>
#include <vtkCellDataToPointData.h>
#include <vtkDataArray.h>
#include <vtkDataSet.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnsignedCharArray.h>
#include <vtkUnstructuredGrid.h>

#include <vtkVisItClipper.h>

#include <avtAccessor.h>
#include <avtIntervalTree.h>
#include <avtMetaData.h>
#include <avtCallback.h>

#include <DebugStream.h>
#include <VisItException.h>
#include <TimingsManager.h>

#include <float.h>

#include <string>
#include <vector>


// ****************************************************************************
//  Method: avtIsovolumeFilter constructor
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Fri Jan 30 14:50:21 PST 2004
//
// ****************************************************************************

avtIsovolumeFilter::avtIsovolumeFilter()
{
}


// ****************************************************************************
//  Method: avtIsovolumeFilter destructor
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Fri Jan 30 14:50:21 PST 2004
//
//  Modifications:
//
// ****************************************************************************

avtIsovolumeFilter::~avtIsovolumeFilter()
{
}


// ****************************************************************************
//  Method:  avtIsovolumeFilter::Create
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Fri Jan 30 14:50:21 PST 2004
//
// ****************************************************************************

avtFilter *
avtIsovolumeFilter::Create()
{
    return new avtIsovolumeFilter();
}


// ****************************************************************************
//  Method:      avtIsovolumeFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: Jeremy Meredith
//  Creation:   January 30, 2004
//
// ****************************************************************************

void
avtIsovolumeFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const IsovolumeAttributes*)a;

    // We need to specify that we want a secondary variable as soon as
    // possible.
    if (strcmp(atts.GetVariable().c_str(), "default") != 0)
    {
        SetActiveVariable(atts.GetVariable().c_str());
    }
}


// ****************************************************************************
//  Method: avtIsovolumeFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtIsovolumeFilter with the given
//      parameters would result in an equivalent avtIsovolumeFilter.
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Fri Jan 30 14:50:21 PST 2004
//
// ****************************************************************************

bool
avtIsovolumeFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(IsovolumeAttributes*)a);
}


// ****************************************************************************
//  Method:  avtIsovolumeFilter::ExecuteSingleClip
//
//  Purpose:
//    Clip against either a low or high scalar.  Convert cell to point
//    data as necessary.
//
//  Arguments:
//    in_ds      the data set to clip
//    val        the value to clip against
//    flip       false if val is a max, true if val is a min
//
//  Programmer:  Jeremy Meredith
//  Creation:    May  6, 2004
//
//  Modifications:
//
//    Hank Childs, Sat Sep 29 11:24:12 PDT 2007
//    Pass in vtkDataArrays to the clipper, not "float *".
//
//    Eric Brugger, Wed Jul 30 18:49:46 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

vtkDataSet *
avtIsovolumeFilter::ExecuteSingleClip(vtkDataSet *in_ds, float val, bool flip)
{
    vtkVisItClipper *clipper = vtkVisItClipper::New();
    clipper->SetInsideOut(flip);

    //
    // Get the scalar array we'll use for clipping; it must be nodal
    //
    vtkCellDataToPointData *cd2pd = NULL;
    if (in_ds->GetPointData()->GetScalars())
    {
        vtkDataArray *s = in_ds->GetPointData()->GetScalars();
        clipper->SetClipScalars(s, val);
    }
    else if (in_ds->GetCellData()->GetScalars())
    {
        //
        // Okay, our active variable was cell-centered.  Recenter it....
        //
        vtkDataSet *temp_ds = (vtkDataSet *) in_ds->NewInstance();
        temp_ds->CopyStructure(in_ds);
        temp_ds->GetCellData()->SetScalars(in_ds->GetCellData()->GetScalars());

        cd2pd = vtkCellDataToPointData::New();
        cd2pd->SetInputData(temp_ds);
        cd2pd->Update();

        vtkDataSet *temporary = cd2pd->GetOutput();

        // Now tell the clipper about it....
        vtkDataArray *s = temporary->GetPointData()->GetScalars();
        clipper->SetClipScalars(s, val);

        // Wait until after the clipping is done to delete 'cd2pd' (which
        // will take 'temporary' with it)
        temp_ds->Delete();
    }
    else
    {
        debug1 << "Could not find any data for isovolume operation\n";
        clipper->Delete();
        EXCEPTION1(VisItException, "No variable was present for the Isovolume");
    }

    //
    // Do the clipping!
    //
    clipper->SetInputData(in_ds);
    clipper->Update();
    vtkDataSet *out_ds = clipper->GetOutput();
    out_ds->Register(NULL);
    clipper->Delete();

    //
    // Free the temporary filter used to convert to point data
    //
    if (cd2pd)
        cd2pd->Delete();

    return out_ds;
}

// ****************************************************************************
// Method: IsovolumeMinMax
//
// Purpose:
//   Determine the min and max for the input data array using an accessor.
//
// Arguments:
//
// Returns:
//
// Note:
//
// Programmer: Brad Whitlock
// Creation:   Sun Apr 22 01:21:24 PDT 2012
//
// Modifications:
//
// ****************************************************************************

template <typename Accessor>
inline void IsovolumeMinMax(double &min, double &max, Accessor access)
{
    min = +FLT_MAX;
    max = -FLT_MAX;
    access.InitTraversal();
    while(access.Iterating())
    {
        double value = access.GetTuple1();
        min = (min < value ? min : value);
        max = (max > value ? max : value);
        ++access;
    }
}

// ****************************************************************************
//  Method: avtIsovolumeFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the Isovolume filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: Jeremy Meredith
//  Creation:   January 30, 2004
//
//  Modifications:
//    Jeremy Meredith, Mon Feb  2 13:13:05 PST 2004
//    Fixed memory leak.
//
//    Jeremy Meredith, Wed May  5 14:56:35 PDT 2004
//    Removed the "3D" from the end of vtkVisItClipper because I made it
//    fully support 2D as well.  I also changed it to only take a single
//    cutoff so the math is more robust, and that required making a "min"
//    pass as well as a "max" pass, so that had to change in this function.
//
//    Jeremy Meredith, Thu May  6 11:38:21 PDT 2004
//    I neglected to create a new scalar array in the event that we do both
//    a min and max pass, so I split most of this routine into a new function
//    and called the new one (ExecuteSingleClip) twice.
//
//    Hank Childs, Wed Oct 20 16:16:15 PDT 2004
//    Basic optimizations -- only do min pass and max pass when the data
//    requires it.  Also pass the dataset through if it is wholly contained
//    within the range and return a NULL dataset if it is outside the range.
//
//    Hank Childs, Wed Nov 17 11:43:53 PST 2004
//    At the end, we try to convert poly-data input back to poly-data output.
//    There was an assumption that after processing the data, it would be
//    in unstructured grid form, which is not true.  Add a check here. ['5640]
//
//    Brad Whitlock, Sun Apr 22 01:21:10 PDT 2012
//    I added double support.
//
//    Gunther H. Weber, Tue May 13 11:51:19 PDT 2014
//    Check whether data set is empty before performing max clip to prevent
//    crash in border case.
//
//    Eric Brugger, Wed Jul 30 18:49:46 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
//    Brad Whitlock, Fri Sep  5 11:09:06 PDT 2014
//    Fix reference counting.
//
//    Eric Brugger, Fri Sep 26 08:49:02 PDT 2014
//    I modified the routine to return a NULL in the case where it previously
//    returned an avtDataRepresentation with a NULL vtkDataSet.
//
//    Eric Brugger, Tue Aug 25 10:13:49 PDT 2015
//    I modified the routine to return NULL if the output data set was NULL.
//
//    Eric Brugger, Mon Feb 13 11:34:09 PST 2023
//    I replaced VTKh with VTKm.
//
// ****************************************************************************

avtDataRepresentation *
avtIsovolumeFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    //
    // Start off by calculating the range of the dataset.
    //
    vtkDataArray *vals = NULL;
    if (in_ds->GetPointData()->GetScalars() != NULL)
        vals = in_ds->GetPointData()->GetScalars();
    else if (in_ds->GetCellData()->GetScalars() != NULL)
        vals = in_ds->GetCellData()->GetScalars();

    if (vals == NULL)
        return in_dr;

    double min, max;
    if(vals->GetDataType() == VTK_FLOAT)
        IsovolumeMinMax(min, max, avtDirectAccessor<float>(vals));
    else if(vals->GetDataType() == VTK_DOUBLE)
        IsovolumeMinMax(min, max, avtDirectAccessor<double>(vals));
    else
        IsovolumeMinMax(min, max, avtTupleAccessor(vals));

    //
    // Check to see if our range is below the min or above the max.  If so,
    // we will have an empty intersection.
    //
    if (max < atts.GetLbound() || min > atts.GetUbound())
    {
        return NULL;
    }

    //
    // Determine if we need to do the min clip or max clip.  Because of the
    // above logic, we can assume that the dataset's max is bigger than
    // the isovolume's lbound and the min is less than the ubound.
    //
    bool doMinClip = false;
    if ((atts.GetLbound() > -1e37) && (min < atts.GetLbound()))
        doMinClip = true;
    bool doMaxClip = false;
    if ((atts.GetUbound() < 1e37) && (max > atts.GetUbound()))
        doMaxClip = true;

    bool doVTKM = VTKmAble(in_dr);
    avtDataRepresentation *out_dr = NULL;
    if (doVTKM && (doMinClip || doMaxClip))
        out_dr = ExecuteData_VTKM(in_dr, {atts.GetLbound(), atts.GetUbound()}, {doMinClip, doMaxClip});
    else
        out_dr = ExecuteData_VTK(in_dr, {atts.GetLbound(), atts.GetUbound()}, {doMinClip, doMaxClip});

    return out_dr;
}


// ****************************************************************************
//  Method: avtIsovolumeFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Indicates the zones no longer correspond to the original problem.
//
//  Programmer: Jeremy Meredith
//  Creation:   February 16, 2004
//
//  Modifications:
//    Kathleen Bonnell, Thu Mar  2 14:26:06 PST 2006
//    Set ZonesSplit.
//
//    Brad Whitlock, Mon Apr  7 15:55:02 PDT 2014
//    Add filter metadata used in export.
//    Work partially supported by DOE Grant SC0007548.
//
// ****************************************************************************

void
avtIsovolumeFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    GetOutput()->GetInfo().GetValidity().ZonesSplit();

    char params[100];
    snprintf(params, 100, "var=%s, min=%lg, max=%lg", atts.GetVariable().c_str(),
        atts.GetLbound(), atts.GetUbound());
    GetOutput()->GetInfo().GetAttributes().AddFilterMetaData("Isovolume", params);
}


// ****************************************************************************
//  Method: avtIsovolumeFilter::ModifyContract
//
//  Purpose:
//      Tell the database that we will need ghost zones.
//
//  Programmer: Hank Childs
//  Creation:   August 11, 2004
//
//  Modifications:
//
//    Hank Childs, Wed Oct 20 17:04:52 PDT 2004
//    Use interval trees to only read in the domains we need.
//
//    Kathleen Bonnell, Mon Aug 14 16:40:30 PDT 2006
//    API change for avtIntervalTree.
//
// ****************************************************************************

avtContract_p
avtIsovolumeFilter::ModifyContract(avtContract_p in_spec)
{
    avtContract_p spec = new avtContract(in_spec);

    const char *varname = NULL;
    if (atts.GetVariable() != "default")
        varname = atts.GetVariable().c_str();
    else
        varname = in_spec->GetDataRequest()->GetVariable();

    //
    // We will need the ghost zones so that we can interpolate along domain
    // boundaries and get no cracks in our isosurface.
    //
    avtDataAttributes &in_atts = GetInput()->GetInfo().GetAttributes();
    bool skipGhost = false;
    if (in_atts.ValidVariable(varname) &&
        in_atts.GetCentering(varname) == AVT_NODECENT)
        skipGhost = true;
    if (!skipGhost)
        spec->GetDataRequest()->SetDesiredGhostDataType(GHOST_ZONE_DATA);

    std::string iso_var(atts.GetVariable());
    if (iso_var == "default")
        iso_var = in_spec->GetDataRequest()->GetVariable();

    avtIntervalTree *it = GetMetaData()->GetDataExtents(iso_var.c_str());
    if (it != NULL)
    {
        double min = atts.GetLbound();
        double max = atts.GetUbound();
        std::vector<int> dl;
        it->GetElementsListFromRange(&min, &max, dl);
        spec->GetDataRequest()->GetRestriction()->RestrictDomains(dl);
    }

    return spec;
}

// ****************************************************************************
//  Method: avtIsovolumeFilter::VTKmAble
//
//  Purpose:
//      Determine if VTKm can be used.
//
//  Programmer: Dave Pugmire
//  Creation:   March 25, 2020
//
//  Modifications:
//    Eric Brugger, Mon Feb 13 11:34:09 PST 2023
//    I replaced VTKh with VTKm.
//
// ****************************************************************************

bool
avtIsovolumeFilter::VTKmAble(avtDataRepresentation *in_dr) const
{
    bool useVTKm = false;
    if (in_dr->GetDataRepType() == DATA_REP_TYPE_VTKM ||
        avtCallback::GetBackendType() == GlobalAttributes::VTKM)
    {
        useVTKm = true;
        vtkDataSet *in_ds = in_dr->GetDataVTK();
        char *var = (activeVariable != NULL ? activeVariable
                                            : pipelineVariable);
        vtkDataArray *pointData = in_ds->GetPointData()->GetArray(var);
        if (pointData == NULL)
        {
            useVTKm = false;
        }
        else if (in_ds->GetDataObjectType() == VTK_RECTILINEAR_GRID)
        {
            vtkRectilinearGrid *rgrid = (vtkRectilinearGrid *) in_ds;

            int dims[3];
            rgrid->GetDimensions(dims);
            if (dims[2] == 1)
                useVTKm = false;
        }
        else if (in_ds->GetDataObjectType() == VTK_STRUCTURED_GRID)
        {
            vtkStructuredGrid *sgrid = (vtkStructuredGrid *) in_ds;
            int dims[3];
            sgrid->GetDimensions(dims);

            if (dims[2] == 1)
                useVTKm = false;
        }
        else if (in_ds->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)
        {
            vtkUnstructuredGrid *ugrid = (vtkUnstructuredGrid *) in_ds;
            vtkIdType nCells = ugrid->GetCells()->GetNumberOfCells();
            vtkUnsignedCharArray *cellTypes = ugrid->GetCellTypesArray();
            unsigned char *ct = cellTypes->GetPointer(0);
            vtkIdType iCell = 0;
            for (; iCell < nCells; ++iCell)
            {
                if (*ct++ != VTK_HEXAHEDRON)
                    break;
            }
            if (iCell != nCells)
                useVTKm = false;
        }
    }

    return useVTKm;
}

// ****************************************************************************
//  Method: avtIsovolumeFilter::ExecuteData_VTK
//
//  Purpose:
//      Perform isoVolume using VTK
//
//  Programmer: Dave Pugmire
//  Creation:   March 25, 2020
//
//  Modifications:
//
// ****************************************************************************

avtDataRepresentation *
avtIsovolumeFilter::ExecuteData_VTK(avtDataRepresentation *in_dr, std::vector<double> bounds, std::vector<bool> clips)
{
    int timerHandle = visitTimer->StartTimer();
    vtkDataSet *in_ds = in_dr->GetDataVTK();
    //
    // Do the clipping!
    //
    vtkDataSet *out_ds = in_ds;
    if(clips[0] && clips[1])
    {
        vtkDataSet *intermediate = NULL;
        intermediate = ExecuteSingleClip(in_ds, atts.GetLbound(), true);
        if(intermediate->GetNumberOfCells() > 0)
        {
            out_ds = ExecuteSingleClip(intermediate, atts.GetUbound(), false);
            intermediate->Delete();
        }
        else
            out_ds = intermediate;
    }
    else if(clips[0])
        out_ds = ExecuteSingleClip(in_ds, atts.GetLbound(), true);
    else if(clips[1])
        out_ds = ExecuteSingleClip(in_ds, atts.GetUbound(), false);
    bool own = out_ds != in_ds;

    //
    // Make sure there's something there
    //
    if (out_ds->GetNumberOfCells() <= 0)
    {
        // We can't delete the out_ds unless we made it and own it.
        // Otherwise, we're deleting the reference out from under the in_ds.
        if(own)
            out_ds->Delete();
        return NULL;
    }

    //
    // If we had poly data input, we want poly data output.  The VTK filter
    // only returns unstructured grids, so convert that now.  Note: we don't
    // necessarily have a ugrid, since it might be that we didn't process the
    // dataset.
    //
    if (in_ds->GetDataObjectType() == VTK_POLY_DATA &&
        out_ds->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)
    {
        vtkUnstructuredGrid *ugrid = (vtkUnstructuredGrid *) out_ds;
        vtkPolyData *out_pd = vtkPolyData::New();
        out_pd->SetPoints(ugrid->GetPoints());
        out_pd->GetPointData()->ShallowCopy(ugrid->GetPointData());
        out_pd->GetCellData()->ShallowCopy(ugrid->GetCellData());
        vtkIdType ncells = ugrid->GetNumberOfCells();
        out_pd->Allocate(ncells);
        for (vtkIdType i = 0 ; i < ncells ; i++)
        {
            int celltype = ugrid->GetCellType(i);
            const vtkIdType *pts;
            vtkIdType npts;
            ugrid->GetCellPoints(i, npts, pts);
            out_pd->InsertNextCell(celltype, npts, pts);
        }
        if(own)
            out_ds->Delete();
        out_ds = out_pd;
        own = true;
    }

    avtDataRepresentation *out_dr = NULL;
    if (out_ds != NULL)
        out_dr = new avtDataRepresentation(out_ds,
            in_dr->GetDomain(), in_dr->GetLabel());

    if (own && out_ds != NULL)
        out_ds->Delete();

    visitTimer->StopTimer(timerHandle, "avtIsovolumeFilter::ExecuteDataTree_VTK");
    return out_dr;
}

// ****************************************************************************
//  Method: avtIsovolumeFilter::ExecuteData_VTKM
//
//  Purpose:
//      Perform isoVolume using VTKm
//
//  Programmer: Dave Pugmire
//  Creation:   March 25, 2020
//
//  Modifications:
//    Eric Brugger, Fri Feb 24 14:57:15 PST 2023
//    I replaced VTKh with VTKm.
//
// ****************************************************************************

avtDataRepresentation *
avtIsovolumeFilter::ExecuteData_VTKM(avtDataRepresentation *in_dr,
                                     std::vector<double> bounds,
                                     std::vector<bool> clips)
{
#ifndef HAVE_LIBVTKM
    return NULL;
#else
    //
    // Get the VTKM data set, the domain number, and the label.
    //
    avtVtkmDataSet *in_ds = in_dr->GetDataVTKm();
    int domain = in_dr->GetDomain();
    std::string label = in_dr->GetLabel();

    if (!in_ds)
    {
        return NULL;
    }

    std::string isoVar(activeVariable != NULL ? activeVariable
                                              : pipelineVariable);

    vtkm::cont::DataSet dataset = in_ds->ds;

    int timerHandle = visitTimer->StartTimer();

    //
    // Recenter the field to the nodes if necessary.
    // FIX_ME: This code is never executed because VTKmAble returns
    // false if the variable is cell associated. This is because the
    // point averaging is flawed. It should create a copy of the field
    // and point average it for doing the isovoluming and then remove
    // the point averaged field afterwards. This way the original cell
    // associated field is left intact.
    //
    bool isCellAssoc = dataset.GetField(isoVar).GetAssociation() ==
                vtkm::cont::Field::Association::Cells;

    if (isCellAssoc)
    {
        vtkm::filter::field_conversion::PointAverage avg;
        avg.SetActiveField(isoVar);
        dataset = avg.Execute(dataset);
    }

    //
    // Execute the clip with field filter.
    //
    vtkm::filter::contour::ClipWithField iso;

    iso.SetActiveField(isoVar);

    //
    // Do lower clip.
    //
    if (clips[0])
    {
        iso.SetClipValue(bounds[0]);
        iso.SetInvertClip(false);
        dataset = iso.Execute(dataset);
    }

    //
    // Do upper clip.
    //
    if (clips[1])
    {
        iso.SetClipValue(bounds[1]);
        iso.SetInvertClip(true);
        dataset = iso.Execute(dataset);
    }

    //
    // Determine if the dataset is empty, and if so, set the output
    // to NULL.
    //
    avtVtkmDataSet *isoOut = NULL;
    if(dataset.GetCellSet().GetNumberOfCells() > 0)
    {
        isoOut = new avtVtkmDataSet;
        isoOut->ds = dataset;
    }

    //
    // Create the output data representation
    //
    avtDataRepresentation *out_dr = new avtDataRepresentation(isoOut, in_dr->GetDomain(), in_dr->GetLabel());

    visitTimer->StopTimer(timerHandle, "avtIsovolumeFilter::ExecuteDataTree_VTKM");

    return out_dr;
#endif
}
